%% process_mesh
function mesh = process_mesh(pde, varargin)
if ~isempty(varargin)
    h = varargin{1};
    if h > (pde.domain(2)-pde.domain(1)) || h > (pde.domain(4)-pde.domain(3)) || h > (pde.domain(6)-pde.domain(5))
        error("h is too large.");
    end
    mesh.Nx = (pde.domain(2)-pde.domain(1))/h;
    mesh.Ny = (pde.domain(4)-pde.domain(3))/h;
    mesh.Nz = (pde.domain(6)-pde.domain(5))/h;
    [mesh.node, mesh.elem] = process_cube(pde.domain, mesh.Nx, mesh.Ny, mesh.Nz);
    [mesh.tri, mesh.bdf] = process_face(mesh);
    [mesh.edge, mesh.elem_ed, mesh.edgeMap] = process_edge(mesh);
    mesh.Np = size(mesh.node,1);
    mesh.Ne = size(mesh.elem,1);
    mesh.Nf = size(mesh.tri,1);
    mesh.Ned = size(mesh.edge,1);
    mesh.Nbd = 6;
else

end
end

%% process_cube
function [Vectex, Tetra] = process_cube(domain, Nx, Ny, Nz)
[X, Y, Z] = ndgrid(linspace(domain(1), domain(2), Nx+1), linspace(domain(3), domain(4), Ny+1), linspace(domain(5), domain(6), Nz+1));
Vectex = [X(:), Y(:), Z(:)];
Npxoy = (Nx+1)*(Ny+1);
idxp = [0 1 Nx+1 Nx+2 Npxoy Npxoy+1 Npxoy+Nx+1 Npxoy+Nx+2];
[Hx, Hy, Hz] = ndgrid(1:Nx, 0:Ny-1, 0:Nz-1);
idx = reshape((Hx + Hy*(Nx+1) + Hz*Npxoy), [], 1);
Hexa = repmat(idxp, size(idx,1), 1)+idx;
Tetra = zeros(6*Nx*Ny*Nz, 4);
t = [5 1 7 8;
    1 5 6 8;
    2 1 6 8;
    1 2 4 8;
    3 1 4 8;
    1 3 7 8];
k = 1;
for i = 1:size(Hexa,1)
    Tetra(k:k+5,:) = reshape(Hexa(i,t), 6, 4);
    k = k+6;
end
% One to five
% Tetra = zeros(5*Nx*Ny*Nz, 4);
% k = 1;
% layer = 1;
% flag_flip = 1;
% for i = 1:size(Hexa,1)
%     row = mod(mod((i-1),Nx),Nx*Ny)+1;
%     col = mod(floor((i-1)/Nx),Ny)+1;
%     if mod(row+col,2) == 0
%         flag = flag_flip;
%     else
%         flag = abs(flag_flip-1);
%     end
%     if flag == 1
%         t = [1 2 3 5;
%             2 3 5 8;
%             3 2 4 8;
%             2 5 6 8;
%             5 3 7 8];
%     else
%         t = [1 5 6 7;
%             4 1 6 7;
%             1 2 4 6;
%             3 1 4 7;
%             6 4 7 8];
%     end
%     Tetra(k:k+4,:) = reshape(Hexa(i,t), 5, 4);
%     k = k+5;
%     if mod(i,Nx*Ny) == 0
%         flag_flip = abs(flag_flip-1);
%         layer = layer+1;
%     end
% end
end

%% process_face
function [Tri, bdf] = process_face(mesh)
elem = mesh.elem;
Ne = size(elem,1);
tri_all = zeros(Ne*4,4);
t = [4 3 2;
    1 3 4;
    4 2 1;
    1 2 3];
k = 1;
for e = 1:Ne
    tri_all(k:k+3,:) = [e*ones(4,1), reshape(elem(e,t), 4, 3)];
    k = k+4;
end
[~, ib, ic] = unique(sort(tri_all(:,2:4),2), "rows", "stable");
k = accumarray(ic,1);
Tri = tri_all(ib,:);
idx = find(k==1);
% Tetra_tri = reshape(ic, 4, [])';
if isfield(mesh, "bdary")
else
    Nhex = mesh.Nx*mesh.Ny*mesh.Nz;
    Hexa_tri = reshape(ic, [], Nhex)';
    tri_face = zeros(6,Nhex);
    t = [7 10;
        8 11;
        14 17;
        4 18;
        13 16;
        2 5];
    k = 1;
    for e = 1:Nhex
        temp = unique(Hexa_tri(e,:), "stable");
        tri_face(:,k:k+1) = reshape(temp(t), 6, 2);
        k = k+2;
    end
    k = ismember(tri_face, idx);
    temp = zeros(size(tri_face));
    temp(k==1) = tri_face(k==1);
    for e = 1:6
        fi = temp(e,temp(e,:)>0);
        Tri(fi,5) = e;
    end
    bdf = Tri(Tri(:,5)>0,:);
end
end

%% process_edge
function [Edge, Tetra_ed, edgeMap] = process_edge(mesh)
elem = mesh.elem;
Ne = size(elem,1);
edge_all = zeros(6*Ne,2);
t = [4 3;
    3 2;
    2 4;
    1 3;
    4 1;
    2 1];
k = 1;
for e = 1:Ne
    edge_all(k:k+5,:) = reshape(elem(e,t), 6, 2);
    k = k+6;
end
[~, ib, ic] = unique(sort(edge_all,2), "rows", "stable");
Edge = edge_all(ib,:);
Tetra_ed = reshape(ic, 6, [])';
E_sort = sort(Edge,2);
edgeMap = containers.Map("KeyType", "char", "ValueType", "any");
for e = 1:size(Edge,1)
    key = sprintf("%d %d", E_sort(e,1), E_sort(e,2));
    edgeMap(key) = e;
end
end